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I. INTRODUCTION 

The recent issue to investigate alternative theories of gravity comes out from Cosmology, Quantum Field Theory and 

Mach's Principle. The initial singularity, the flatness and horizon problems [1] point out that Standard Cosmological 

Model [2], based on General Relativity (GR) and Particle Standard Model, fails in describing the Universe at extreme 

regimes. Besides, GR does not work as a fundamental theory capable of giving a quantum description of spacetime. 

&Q. Due to these reasons and to the lack of a definitive Quantum Gravity theory, alternative theories of gravitation have 

been pursued in order to attempt, at least, a semi-classical approach to quantization. In particular, Extended Theories 

r*-\ , of Gravity (ETGs) face the problem of gravitational interaction correcting and enlarging the Einstein theory. 

^ ' The general paradigm consists in adding, into the effective action, physically motivated higher-order curvature 

C*") | invariants and non- minimally coupled scalar fields [3, 4]. 

^^•O The interest of such an approach in early epoch cosmology is due to the fact that ETGs can "naturally" reproduce 

inflationary behaviors able to overcome the shortcomings of the Standard Cosmological Model and seems also capable 
of matching with several observations. 

From another viewpoint, the Mach Principle gives further motivations to modify GR stating that the local inertial 
frame is determined by the average motion of distant astronomical objects [5]. As a consequence, the gravitational 
coupling can be scale-dependent. This means that the concept of inertia and the Equivalence Principle have to be 
revised since there is no a priori reason to restrict the gravitational Lagrangian to a linear function of the Ricci scalar 
R, minimally coupled with matter [6—11]. 

Very recently, ETGs are playing an interesting role to describe today's observed Universe. In fact, the impressive 
amount of good quality data of last decade seems to shed new light into the effective picture of the Universe. Type 
la Supernovae (SNela) [12], anisotropics in the CMBR [13], and matter power spectrum derived from wide and deep 
galaxy surveys [14] represent the strongest evidences for a radical revision of the Cosmological Standard Model also 
at recent epochs. 

Specifically, the Concordance A COM Model is showing that baryons contribute only for ~ 4% to the total matter - 
energy budget, while the cold dark matter (CDM) represents the bulk of the clustered large scale structures (~ 25%) 
and the cosmological constant A plays the role of the so called "dark energy" (~ 70%) [15]. 

Although being the best fit to a wide range of data [16], the ACDM model is affected by strong theoretical 
shortcomings [17] that have motivated the search for alternative models [18, 19]. 

Dark energy models mainly rely on the implicit assumption that Einstein's GR is the correct theory of gravity 
indeed. Nevertheless, its validity on large astrophysical and cosmological scales has never been tested but only 
assumed [20], and it is therefore conceivable that both cosmic speed up and missing matter are nothing else but 
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signals of a breakdown of GR. In this sense, GR could fail in giving self-consistent pictures both at ultraviolet scales 
(early universe) and at infrared scales (late universe). 

Following this line of thinking, the "minimal" choice could be to take into account generic functions f(R) of the 
Ricci scalar R. However, such an approach can be encompassed in the ETGs being the minimal extension of GR. 
The task for this extended theories should be to match the data under the "economic" requirement that no exotic 
dark ingredients have to be added, unless these are going to be found with fundamental experiments [21]. This is the 
underlying philosophy of what are referred to as f(R)-gravity (see [19, 22, 23] and references therein). 

Although higher order gravity theories have received much attention in cosmology, since they are naturally able to 
give rise to the accelerating expansion (both in the late and in the early universe [24]), it is possible to demonstrate 
that f(R) theories can also play a major role at astrophysical scales. In fact, modifying the gravity Lagrangian affects 
the gravitational potential in the low energy limit. Provided that the modified potential reduces to the Newtonian 
one on the Solar System scale, this implication could represent an intriguing opportunity rather than a shortcoming 
for f(R) theories. In fact, a corrected gravitational potential could offer the possibility to fit galaxy rotation curves 
without the need of huge amounts of dark matter [25-31]. In addition, it is possible to work out a formal analogy 
between the corrections to the Newtonian potential and the usually adopted galaxy halo models which allow to 
reproduce dynamics and observation without dark matter [27]. 

However, extending the gravitational Lagrangian could give rise to several problems. These theories could have 
instabilities [32, 33] and ghost -like behaviors [34-36], and they have to be matched with the low energy limit exper- 
iments which fairly test GR. Besides, these theories should also be compatible with early universe tests such as the 
formation of CMBR anisotropics, Big Bang Nucleosynthesis [37], and Baryogenesis [38, 39]. 

Actually, the debate concerning the weak field limit of /(i?)-gravity is far to be definitive. In the last few years, sev- 
eral authors have dealt with this matter with contrasting conclusions, in particular with respect to the Parameterized 
Post Newtonian (PPN) limit [40, 42]. 

In summary, it seems that the paradigm to adopt /(i?)-gravity leads to interesting results at cosmological, galactic 
and Solar System scales but, up to now, no definite physical criterion has been found to select the final f(R) theory 
(or class of theories) capable of matching the data at all scales. Interesting results have been achieved in this line of 
thinking [21, 43-46] but the approaches are all phcnomenological and are not based on some fundamental principle 
as the conservation or the invariancc of some quantity or some intrinsic symmetry of the theory. Furthermore, as it 
was shown in [32] , in alternative theories of gravity, it is important to understand the background before exploring 
other bounds, such as anisotropics in the CMBR. For this goal it is essential to try to find exact analytical solutions 
for the f(R) theories, and, only after this, study more in detail the possible evolutions compatible with our data (e.g. 
solar system and CMBR bounds). 

In some sense, the situation is similar to that of dark matter: we know very well its effect at large astrophysical 
scales but no final evidence of its existence has been found, up to now, at fundamental level. In the case of f(R)- 
gravity, we know that the paradigm is working: in principle, the missing matter and accelerated cosmic behavior can 
be addressed taking into account gravity (in some extended version), baryons and radiation but we do not know a 
specific criterion to select the final, comprehensive theory. 

In this paper, we want to address the following issues: i) Is there some general principle capable of selecting 
physically motivated f(R) models? ii) Can conserved quantities or symmetries be found in relation to specific f(R) 
theories? Hi) Can such quantities, if existing, give rise to viable cosmological models? 

In this paper, following the so called Noether Symmetry Approach (see [7, 47, 48], we want to seek for viable f{R) 
cosmological models. As we will see, the method is twofold: from one side, the existence of symmetries allows to solve 
exactly the dynamics; from the other side, the Noether charge can always be related to some observable quantity. 

The layout of the paper is the following. In Sec. II, we sketch the dynamics of f(R) gravity in the metric approach 
and derive the Friedmann-Lemaitre-Robertson- Walker (FLRW) cosmological equations. Sec. Ill is devoted to the 
general discussion of the Noether Symmetry Approach by which it is possible to find out conserved quantities and 
then symmetries which allow to exactly solve a dynamical system. In Sec. IV, we apply the method to the f(R) 
cosmology. In Sec.V, we give a detailed summary of the exact solutions discussing them in presence or in absence of 
the Noether charge. Sec. VI is devoted to the discussion and the conclusions. 

II. f(R) GRAVITY AND COSMOLOGY 

The action 

S= [ d A x^^f(R) + S m , (1) 



describes a theory of gravity where f(R) is a generic function of the Ricci scalar R. GR is recovered in the particular 
case f(R) = —R/l6irG, and S m is the action for a perfect fluid minimally coupled with gravity . 

This action, in general, leads to 4th order differential equations for the metric since the field equations are 

Jr. R^v - 2 i 9\a> ~ fR;nv + g^v Of R — — 2 T™ , (2) 

where a subscript R denotes differentiation with respect to R and T™ is the matter fluid stress-energy tensor. 
Defining a curvature stress - energy tensor as 

T%r = -^ j^„ [f(R) - Rf R (R)} + f R {Ry a(i {9^gp v - g^gap)} , (3) 



Eqs.(2) can be recast in the Einstein - like form: 



G„ v = R^ - - 9lw R = T$T + T™/f R {R) (4) 



where matter non - minimally couples to geometry through the term 1/ f R (R). It is known that these theories can be 
mapped to a scalar-tensor theory. However, there are two points which should be noticed. First, the two theories 
might have different quantum descriptions, as they only coincide on the classical solutions. Furthermore, the two 
theories are classically equivalent if the Brans-Dicke parameter (wbd) exactly vanishes and if the scalar field possesses 
a suitable potential. This fact is related to the second point: in the literature, the Brans-Dickc field is commonly 
taken as a light scalar field for which the local gravity constraint fixes the Brans-Dickc parameter to be greater than 
40000. This bound is usually considered when studying Brans-Dicke theories. However, for the f(R) theories, since 
wbd = 0, this is not the case, and the presence of a non-negligible potential is essential in order to give an explicit 
mass to the gravitational scalar degree of freedom. Once one has the solution H(t) (and consequently R(t)) for a 
given f(R), the scalar field is defined as $(£) = -f R (t), and its potential is £/($(£)) = R(t) f R (t) - f(R(t)). An 
example showing this link between scalar-tensor theories and f(R) gravity is given in the appendix for one solution 
which will be found explicitly later on. 

In order to derive the cosmological equations in a FLRW metric, one can define a canonical Lagrangian C = 
C(a,a,R,R), where Q = {a,R} is the configuration space and TQ = {a,a,R,R} is the related tangent bundle 
on which C is defined. The variable a(t) and R(t) are the scale factor and the Ricci scalar in the FLRW metric, 
respectively. One can use the method of the Lagrange multipliers to set R as a constraint of the dynamics. Selecting 
the suitable Lagrange multiplier and integrating by parts, the Lagrangian L becomes canonical. In our case, we have 



5 = 2tt 2 / dt<f{ f(R)- X 
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where a is the scale factor scaled with respect to today's value (so that a = a/ao and a(to) = 1); p m o and p r o represent 
the standard amounts of dust and radiation fluids as, for example, measured today; finally k — fc/ag, where k — 0, ±1. 
This choice for a, makes it dimensionless, and it also implies that [k] = [R] = M 2 , whereas [/] = [p r o] = M A - It is 
straightforward to show that, for f(R) = —R/WttG — paq, one obtains the usual Fricdmann equations. 

The variation with respect to R of the action gives A = f R . Therefore the previous action can be rewritten as 



S = 2tt 2 / dta 3 if-f R 
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and then, integrating by parts, the point-like FLRW Lagrangian is 

C = a 3 (/ - f R R) + 6 a 2 f RR R a + 6 f R a a 2 - 6k f R a - p m0 - p r0 /a , (7) 

which is a canonical function of two coupled fields, R and a, both depending on time t. The total energy Ec, 
corresponding to the 0, 0-Einstcin equation, is 

E c = 6 f RR a 2 aR + 6 f R aa 2 - a 3 (f - f R R) + 6k f R a + p m0 + ^ =0. (8) 

a 
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As we shall see later, it is convenient to look for parametric solutions in the form \H(a), /(i?(a))] , so that f R = f'/R', 
where a prime denotes differentiation with respect to the time-parameter a. We also have that, if R ^ constant, 
Irr R = dfn/dt = a H f' R = a H [f'/R 1 — f R"/R' ], so that the Friedmann equation can be rewritten as 
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The equations of motion for a and R are respectively 
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where H = a/a is the Hubble parameter. Considering R and a as the variables, we have, for consistency (excluding 
the case f RR — 0), that R coincides with the definition of the Ricci scalar in the FLRW metric. Geometrically, this is 
the Euler constraint of the dynamics. Using (10), only one of the equations (8), and (11) is independent because of 
the Bianchi identities, as these equations correspond to the first and second modified Einstein equations, and matter 
is conserved. Equivalently, after multiplying equation (11) by a 2 a, and using (10), one can integrate (11) to find (8). 
Furthermore, as we will show below, constraints on the form of the function f(R) and, consequently, solutions of the 
system (8), (10) can be achieved by asking for the existence of Noether symmetries. Such solutions will also solve 
equation (11) automatically. On the other hand, the existence of the Noether symmetries guarantees the reduction 
of dynamics and the eventual solvability of the system. 



III. THE NOETHER SYMMETRY APPROACH 



Solutions for the dynamics given by (7) can be achieved by selecting cyclic variables related to some Noether 
symmetry. In principle, this approach allows to select /(-R)-gravity models compatible with the symmetry so it can 
be seen as a physical criterion since the conserved quantities are a sort of Noether charges. Therefore such a criterion 
might be to look for those f(R) which have cosmological Noether charge. Although this criterion somehow "breaks" 
Lorentz-invariance because we need the FLRW background to formulate it, however Lorentz-invariance is evidently 
broken in our universe by the presence of the CBMR radiation which, by itself, fixes a preferred reference frame. 

In general, the Noether Theorem states that conserved quantities are related to the existence of cyclic variables 
into dynamics [49-51]. 

Let C(q l ,q l ) be a canonical, non-degenerate point-like Lagrangian where 
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d 2 C 
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(12) 



with Hij the Hessian matrix related to C and a dot denotes differentiation with respect to the affine parameter A. 
The dot indicates derivatives with respect to the affine parameter A which, in our case, corresponds to the cosmic 
time t. In analytical mechanics, C is of the form 



£ = T(q,q)-V(q) 



(13) 



where T and V are the "kinetic" and "potential energy" respectively. T is a positive definite quadratic form in q. 
The energy function associated with C is 



E r 






(14) 



which is the total energy T + V . In any case, E^ is a constant of motion. Since our cosmological problem has a finite 
number of degrees of freedom, we are going to consider only point-transformations. Any invcrtiblc transformation of 
the "generalized positions" Q l = Q l (q) induces a transformation of the "generalized velocities" such that 



m = %<? 



(15) 



the matrix J = ||9(5 4 /^9"'ll i s the Jacobian of the transformation on the positions, and it is assumed to be nonzero. 
The Jacobian J of the induced transformation is easily derived and J 7^ — >• J 7^ 0. In general, this condition is not 
satisfied in the whole space but only in the neighbor of a point. It is a local transformation. 

A point transformation Q % — Q l (q) can depend on one (or more than one) parameter. As starting point, we 
can assume that a point transformation depends on a parameter e, i.e. Q % — Q l (q, e), and that it gives rise to a 
one-parameter Lie group. For infinitesimal values of e, the transformation is then generated by a vector field: for 
instance, d/dx is a translation along the x axis, x(d/dy) — y(d/dx) is a rotation around the z axis and so on. The 
induced transformation (15) is then represented by 

x -^w + {pto)w- (16) 

X is called the "complete lift" of X [51]. A function F(q, q) is invariant under the transformation X if 

^««'<">g + (sr' (q) )S- ' (17) 

where L-%_F is the Lie derivative of F. Specifically, if Ly^C = 0, X is a symmetry for the dynamics derived by C 
As we shall see later on, we will look for a sufficient condition on the form of f(R) in our Lagrangian, which allows 
Lx£ = to vanish. 

Let us consider now a Lagrangian C and its Euler-Lagrange equations 

d dC dC _ 

Let us consider also the vector field (16). Contracting (18) with the a l 's gives 

0. (19) 



Being 



from (19), we obtain 
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The immediate consequence is the Noether Theorem which states: 
If Lx£ = 0, then the function 

So = 4£, (22) 

is a constant of motion. 

Some comments are necessary at this point. Eq.(22) can be expressed independently of coordinates as a contraction 
of X by a Cartan one-form 

dc = ^ . (23) 

For a generic vector field Y = y l d/dx l , and one-form j3 = fiidx 1 , we have, by definition, iy(3 = y % j3i- Thus Eq.(22) 
can be written as 

irfc = So • (24) 

By a point-transformation, the vector field X becomes 

^'^^HsHi' (25) 



We see that X' is still the lift of a vector field defined on the "space of positions." If X is a symmetry and we choose 
a point transformation such that 

ixdQ 1 = 1 ; i x dQ* = i^l, (26) 

we get 

Thus Q 1 is a cyclic coordinate and the dynamics results reduced [49, 50]. 

Furthermore, the change of coordinates given by (26) is not unique and then a clever choice could be very important. 
In general, the solution of Eq.(26) is not defined on the whole space. It is local in the sense explained above. Besides, 
it is possible that more than one X is found, e.g. Xi, X2. If they commute, i.e. [Xi, X2] = 0, then it is possible to 
obtain two cyclic coordinates by solving the system 

ix 1 dQ 1 = 1; ix 2 dQ 2 = 1; i^dQ 1 = 0; % ^ 1; i^dQ 1 = 0; % ^ 2 . (28) 

The transformed fields will be d/dQ 1 , d/dQ 2 . If they do not commute, this procedure is clearly not applicable, since 
commutation relations are preserved by diffcomorphisms. If the relation X3 = [Xi, X2] holds, also X3 is a symmetry, 
being L-^ a C = L^L^C — Lx 2 Lx 1 £ = 0. If X 3 is independent of Xi, X 2 , we can go on until the vector fields close 
the Lie algebra. The usual approach to this situation is to make a Legendre transformation, going to the Hamiltonian 
formalism, and then derive the Lie algebra of Poisson brackets. 

If we seek for a reduction of dynamics by cyclic coordinates, the procedure is possible in the following way: i) we 
arbitrarily choose one of the symmetries, or a linear combination of them, searching for new coordinates where, as 
sketched above, the cyclic variables appear. After the reduction, we get a new Lagrangian A(Q); ii) we search again 
for symmetries in this new configuration space, make a new reduction and so on until possible; in) if the search fails, 
we try again by another of the existing symmetries. 

Let us now assume that C is of the form (13). As X is of the form (16), L^C will be a homogeneous polynomial of 
second degree in the velocities plus a inhomogeneous term in the q % . Since such a polynomial has to be identically zero, 
each coefficient must be independently zero. If n is the dimension of the configuration space, we get {1 + n(n + l)/2} 
partial differential equations. The system is overdetermined, therefore, if any solution exists, it will be expressed in 
terms of integration constants instead of boundary conditions. It is also obvious that an overall constant factor in the 
Lie vector X is irrelevant. In other words, the Noether Symmetry Approach can be used to select functions which 
assign the models and such functions (and then the models) can be physically relevant. 

Considering the specific case which we are going to discuss, the f(R) cosmology, the situation is the following. The 
configuration space is Q = {a,R} while the tangent space for the related tangent bundle is TQ = {a,a,R,R}. The 
Lagrangian is an application 

C : TQ — ► n (29) 

where !K is the set of real numbers. The generator of symmetry is 

d d d ■ d 

da oR da dR 

As discussed above, a symmetry exists if the equation L-^C = has solutions. Then there will be a constant of motion 
on shell, i.e. for the solutions of the Eulcr equations, as stated above equation (22). In other words, a symmetry exists 
if at least one of the functions a or f3 in Eq.(30) is different from zero. As a byproduct, the form of f(R), not specified 
in the point-like Lagrangian (7), is determined in correspondence to such a symmetry. 

IV. NOETHER SYMMETRIES IN f(R) COSMOLOGY 

For the existence of a symmetry, we can write the following system of equations (linear in a and /3), 

f R {a + 2ad a a) + a f RR (f3 + ad a l3) = (31) 

a 2 f RR d R a = (32) 

2f R d R a + f RR (2a + ad a a + ad R (3)+af3f RRR = 0, (33) 



obtained setting to zero the coefficients of the terms a 2 , R 2 and aR in Ly^C = 0. In order to make L-^C 
we will also look for those particular /'s which, given the Eulcr dynamics, also satisfy the constraint 



3a (/ - R.fn) - a/3R f RR - (af R + a /3 f RR ) 
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(34) 



This procedure is different from the usual Noether symmetry approach, in the sense that now L-^C = will be solved 
not for all dynamics (which solve the Euler-Lagrange equations), but only for those / which allows Euler solutions 
to solve also the constraint (34). Imposing such a constraint on the form of / will turn out to be, as we will show, 
a sufficient condition to find solutions of the Euler-Lagrange equation which also possess a constant of motion, i.e. a 
Noether symmetry. As we shall see later on, the system (31), (32) and (33) can be solved exactly. Having a non-trivial 
solution for a and f3 for this system, one finds a constant of motion if also the constraint (34) is satisfied. In fact, with 
these a and /?, only those Euler-Lagrange solutions which also satisfy equation (34) will have a constant of motion. 
However, this will not happen for all f(R)'s. The task will be to find such forms of /. 

A solution of (31), (32) and (33) exists if explicit forms of a, j3 are found. If, at least one of them is different from 
zero, a Noether symmetry exists. 

If fRR 7^ 0, Eq.(32) can be immediately solved being 

a = a{a) . (35) 

The case f RR — is trivial since corresponds to the standard GR. We can rewrite Eqs.(31) and (33) as follows 



./'„> ( a + 2a ^ J + a f RR (fi + a d a p) = 
Irr ( 2 a + a — + a d R /3 J + a j3 f RRR = . 



(36) 
(37) 



Since / = f(R), then df /da = 0, even in the case we consider R = R(a), it is possible to solve equation (37), by 
writing it as 



0j<\ '>fnn) = -//;/; (2- + ^ 
a da 



whose general solution can be written as 
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Therefore one finds that Eq. (36) gives 
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which has solution 
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(38) 



(39) 



(40) 



a = c\ a -\ and h = — 

a a 



(41) 



where, being a dimensionless, c\ and C2 have the same dimensions. We can further fix a to be dimensionless, this 
fixes the dimensions of j3 to be [0] = M 2 . Then also [c] = M 2 , so that we have 
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Irr a f RR 
We can now use the expressions for a and j3 into Eq.(34) as follows 

3 a (a a 2 + c 2 ) / - c (a 2 R + 6k) (cia 2 + c 2 ) p r0 



Ir — 



2a(c2-R — 6cik) 



2a 4 (c 2 i?-6ciK) ' 



(42) 



(43) 



if C2 R — 6 k c\ 7^ 0. It is clear that, for a general /, it will not be possible to solve at the same time the Euler-Lagrange 
equation and this constraint. Therefore we have to use the Noether constraint in order to find the subset of those 



/ which make this possible. As we shall see later, it is convenient to look for a parametric solution in the form 
[H(a), /(i?(a))] . In this case, since f R = f'/R', the Noether condition corresponds to the following ODE 

f'(a) _ 3 a (ci a 2 + e 2 ) /(a) - c (a 2 R(a) + 6k) (c lQ 2 + c 2 ) p r0 

R'(a)~ 2a(c 2 R(a) - 6cm) + 2a 4 (c 2 i?(a) - 6ci«) ' 

It should be noted that this change of variable is defined only if R' ^ 0, that is if R is not constant during the evolution. 
When this happens Eq. (34) or (45) sets a = ciq = constant, which corresponds to an uninteresting solution. 

Any Euler-Lagrange solution, by definition, satisfies the Einstein equations. However we will show that there are 
forms of f(R), for which a subset of those solution will also be a Noether solution. In fact, Eq.(43) can also be 
rewritten as 

Cl a 2 {p r0 + 3 a 4 / + 12 k a 2 f R ) + c 2 [p r0 + a 4 (3 / - 2 R f R )] = ca 3 (a 2 R + 6k) . (45) 

Therefore we look for a family of solutions that, being a Noether symmetry, gives a class of f(R) models. 
This symmetry implies the existence of the following constant of motion 

Oi (6 f RR a 2 R + 12 f R a a) + j3 (6 f RR a 2 a) = 6 /ijj = constant, (46) 

where po has the dimensions of a mass. Equation (46) can be recast in the form 

,, - jrr R - —, a— r H 3— — JR H 5" — H ' i 47 ) 

at a(cia z +C2) c\ a z + c 2 C\ a A + c 2 

or, using the time-parameter a 

nH (rA (HM iM^im (" 2 ci-c 2 )H(a)f'(a) ^ ca 

[ >\R'(a) i?'(a) 2 J (cia 2 + c 2 )i?'(a) a( Cl a 2 + c 2 ) Cl a 2 + c 2 lJ ' lJ 

Once Eq. (44) is solved, because the Noether constraint is satisfied, the solution \H(a), f(R(a))] will automatically 
solve also (48) for a particular ^o- Equation (46) can be used to reduce the order of the Friedmann equation. In fact, 
writing Eq.(8) as 

/ -6f RR RH-6f R H 2 -f R (R+^l)-P^-P^=0, (49) 



we have 



12 ci a 2 2 / 6k\ 6co_ 2 6^gg fao Pro 
—. jRti - } R \R + — M 5— ii = — — r H 3- + — j- 

ci a z + C2 \ a z / ci a^ + C2 0(018^ + 02) a 6 a 4 



where /# is given by (43). We will use this relation in order to find out exact cosmological solutions. Namely, we will 
search for solutions depending on the constant of motion /zo determined by the Noether symmetry. 

V. EXACT COSMOLOGICAL SOLUTIONS 

In order to find out exact cosmological solutions, let us discuss the Noether condition Eq.(45) and the dynamical 
system (8), (10) with respect to the values of the integration constants Ci i2 , the structural parameters k, p r oi Pmo and 
the Noether charge p,Q. Beside cosmological solutions, also the explicit form of f(R) will result fixed in the various 
cases. As we shall see later on, analytical solutions can be easily found for the case when both c and po vanish at the 
same time. Therefore in all this section, except one subsection, we will set c = 0. 

A. Case C\ — 

In this case, the Noether condition (45) reduces to 

2R.f R -3.f=^-. (51) 



1. Vacuum and pure dust case 



In vacuum, or in the presence of dust only, i.e. p r o = 0, we find 



/ - /o ( -5- 



3/2 



(52) 



This solution, for the vacuum case p r o — p m o = 0, has been already found [48]. The absence of a ghost imposes that 
fji < 0, i.e. /o > since i?o < 0. In the case of dust and no radiation (p m o =/= 0,p r o = 0), one can substitute Eq.(52) 
into (50) to find 



R_ 

Ro 



3/2 



18k 



R_ 

Rq 



1/2 



12pf t H _ 2p ma 
c 2 a f a a 3 f 



(53) 



1. k = 0. In this case, for consistency, we need the right hand side of (53) to be positive. If p = (case for which 
analytical solutions could be given), this is impossible as /o > 0, therefore there is no ghost-free solution. For 
the more general case /Xq/c 2 < 0, there could be a physical solution: the non-linearity of the equations does not 
allow us to find analytical solutions for this case. Nevertheless, solutions (to be found numerically) may still 
exist. 

2. k =/= 0. The Ricci scalar can be found as the solution of Eq. (53). For po — 0, we have a cubic equation in 
{R/Rq) 1 / 2 , for which a real solution always exists (which may not be positive though). Looking at equation 
(53), the case /Uo = 0, k = — 1 has no ghost-free solutions (/o < 0). Also the case ^0 = 0, k = 1 has no solution, 
because we have 



(54) 



where we have defined the constant 




Bo = V/oVM + 216 /0 « 3 Ro ~ foPmoRl 



(55) 



which implies that (fo/pmo) 2 (k/Ro) 3 > —1/216. If so, then, since Rq < 0, B n > 0. However, this would lead 
to a negative value for {R/Rq) 1 / 2 . 



2. Dust and radiation case 



In this case we have 



f R 



3 / 



PrO 



2 R 2a 4 R' 
Once again, in order to have fa < 0, and R < during the evolution of the universe one requires 



f >- — 



PrO 

3 a 4 



If we substitute the expression for fn into the reduced Friedmann Eq.(50) we find 

12 /Jq a HR 6k p r Q ipr-oR 2pmoR 



/ = - 



c 2 (R a 2 + 18k) a 4 (Ra 2 + 18k) a 2 (Ra 2 + 18k) a (Ra 2 + 18k) 



(56) 



(57) 



(58) 



This relation gives / as a function of a being R — R(a). It has to be c 2 ^ otherwise the Noether condition becomes 
trivial. This expression can be inserted back into (56). Assuming R = R(a) as a monotonic function of a, one finds 
that fn = (df/da)/(dR/da), and equation (51) becomes a differential equation for R(a), which can be written as 



Rl 



(i 



x {-R 2 l2a 3 (H-aH')ii 3 ] 



a 3 (18a 3 Hpl + 4c 2 /5 r0 + 3ac 2 p m o) (Ra 2 + 6k) 
+ c 2 (2p r0 + ap m0 )]a 4 + 6KR[6a 3 pl (H + aH') - c 2 (4/9 r0 + ap m0 )]a 2 - 72c 2 K 2 p r0 } . 



(59) 



10 

where the prime denotes differentiation with respect to the scale factor a. Eq.(59) can be further rewritten as a second 
order differential equation in H(a), by using equation (10), 



K 



R= -12 H 2 -QaHH' -6 — . (60) 

a 



Substituting (60) into (59) one finds 

H" = , rr , /1fl , rr o - r x {24aK 2 p? l + H[a 2 {a 2 (6a 3 Hp 3 + 4c 2 p r o 

a 4 iJ 2 (I8a 3 ff^ + 4c 2 p r0 + 3ac 2/ 9 m0 ) l ° L l V ° 

+ 3ac 2 p m0 )H' 2 + a[12anpl + H(78a 3 Hp 3 + 32c 2 p r0 + 21ac 2 p m o)]H' 

+ 12H[2anp 3 + H(2a 3 Hp 3 + 2c 2/ o r0 + ac 2 p m0 )}} - 8c 2 K Pr0 ]} . (61) 

This differential equation selects those f(R) models which satisfy, at the same time, both the Friedmann equation and 
the Noether condition. It has to be stressed that, having chosen a as the time variable, finding the iJ(a)'s which solve 
(61) uniquely fixes the metric tensor. Hence, H (a) represents a fully solved exact solution for the Einstein equations. 
Of course, if one wants to know the link between a and the proper time, a = a(t), one needs to find the integral 
t = jda/{aH). 

The case po = is interesting as it allow us to find analytical solutions, as the differential equation be- 
comes (2nd order and) linear for the variable H 2 . In this case, the solution of the equation will be a family 
H = H{a, di,d 2l c 2l po, n, p r o, p m o), where di ]2 are two constants coming from the integration of Eq.(61). In turn, by 
using Eq.(60), it is possible to define a function R = R(a, d\, d 2 ,c 2 , po, Pro, Pmo)> which can then be substituted into 
Eq.(58) in order to find the explicit parametric form of f(R), i.e. / = /(a, d\ 1 d 2 , c 2 , po, p r o, p m a)- In other words, we 
find the explicit parametric form for f(R) where the parameter used to describe the f(R) is the scale factor a (see 
also [21] for a comparison with observations. However, in that case, the adopted f(R) models were constructed by 
phenomenological considerations and not derived from some first principle, as the existence of symmetries as discussed 
here). 

We can distinguish some relevant cases. 

1. k = 0, po = 0. In this case, by exactly integrating equation (61), we find 

rj 2 , d x + 8 a p r0 + 3 p m0 a 2 

H =d 2 , (62) 

a* 

where d\^ 2 are integration constants, with [di] — M A and [d 2 ] = M~ 2 . This expression for H(a) together with 
(58) and (60) form a solution for the set of ODE's (9), and (44), so that Eq. (48) is satisfied giving po = 0. 
Although this solution is analytical it cannot be accepted because it allows for a negative Newton constant. In 
fact, equation (57) cannot be satisfied by equation (58) if k = 0, po = 0. However the non-linear case po/c 2 < 
could actually lead to physical solutions (to be discussed elsewhere in a forthcoming paper). For the same 
reason, also the case k = — 1, po = should be rejected. 

2. k = 1, po = 0. As far as R < —l&n/a 2 , the second term in the l.h.s. of equation (58) becomes positive, allowing 
for the possibility of finding a physical solution. The integration of (61) leads to 

H 2 = (V2 dl - 3 44^) 1 + (Wo - ^) -4 + 3Jl ^ , (63) 



9Pom / a V 3p, 

with [di] = M 2 , and [d 2 ] = M~ 2 . In order to find d\ and d 2 one can fit this formula with the standard Friedmann 
equation of GR with only matter, radiation and curvature. Therefore, one has to consider 



9pQm 



V2 dl -°-^^ = H 2 nf , (64) 

16/CVo K 



8d 2 p r0 -^^ = H 2 n^ , (65) 

3d 2 pm0 - Ho^to, (66) 



but this system admits no solutions as one finds 



2 ^0 ^fcO " 16 ~^~ H ^m0 < ° (67) 

PrO 



using today's data [52]. 
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B. Case C2 = 

In this case, the Noether condition (45) reduces to 

p r0 + 3a 4 f + 12 K a 2 f R = 0. (68) 

1. Vacuum and dust only case 

In this case we have p r o = 0, and a flat universe cannot be solution as one would obtain / = 0. Considering k =/= 
one finds 

Since Jr<0 then / is positive when k < and viceversa. Substituting this into the Friedmann equation one finds 

{a 3 ci [(12H 2 + R)a 2 + 10k}} f = 4k (6Hpl + ciPmo) • (70) 

Restricting ourselves only to the study of the simple and linear case of a vanishing p , we can distinguish two cases 
1- PmO = 0, Ho = 0. In this case one needs to impose 

R = -12H 2 -10-^, (71) 

ar 

which, together with the definition of R, gives 

H 2 = 2d 1 -^, (72) 

where d\ is a constant of integration with dimensions M 2 . This behavior describes a universe with only a 
cosmological constant and curvature. Equation (68) can now be solved for /(a) giving 



t d2 A 

J = — = «2 

a 



R + 24d x 



1/2 

(73) 



2 k 

where d 2 is a constant of integration with dimensions M . 
2- Pmo 7^ 0,/io = 0. In this case the Friedmann equation and (69) give 

/ = ^2 (74) 

J (12H 2 + R) a 5 + IOkci 3 K ' 

Substituting this expression in (69), and using the definition for R in terms of H(a) one finds a linear 2nd order 
differential equation in H 2 (a), which has solution 



where d\^ are integration constants, and [d\] = [d 2 ] = M 2 . Therefore one has 



R = -24d 2 - Hj , (76) 

f = ~ E ?- (77) 

6 adi 
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2. Radiation and dust case 



Also in this case, we have three possibilities, according to the values of k. 
1. k — 0. In this case one finds that 



Therefore we have 



'-&■ < 78) 

• fe -^-3^R7- ( 79 ) 

A well-behaved background evolution requires, with our conventions, R' > 0, so that /r > 0. This means a 
negative effective Newton constant, i.e. the solution cannot be accepted. 

2. k 7^ 0. In this case, using equation (68) one finds 

, _ Pro fa 2 

fR -~T2^-^> (80) 

and then using Friedmann equation (50) one can solve for /, as follows 

-Ci (I2g 2 + R) p r0 a 2 + 12k (QHpl + Cip^o) a + &cinp r0 
■' ~ 3a 4 ci [(12H 2 + R)a 2 + 10k] ' ( j 

By plugging this relation into the Noether condition (68), and using the definition of R in terms of H,H', and 
a, one finds the following differential equation for H(a) 



aH 
H" = 



— (l8aH/iQ + 3acip TO o + 4ci/O r o) H' a 4 — 3 (aH (30aHpo + 5acip m o + 8c\p r o) — 4k/Xq) H' a 2 



a 5 H 2 (18aHp,Q + 3acip m Q + 4ci/9 r o) 
4k (SaHpl + acipmo + 2cip r0 )] - 8k 2 pi 

+ . Sr ,nio-rr..3 I o , T. I N ' i 82 ) 



a 5 _ff 2 (18aHpQ + 3acip m o + 4c\p, 



rO) 

In the case po = 0, p m o ^ 0, this differential equation can be exactly integrated to give 

2 _ 256k/^ _ 16k P 2 8dip ra _ 2k_ 3p m0 rfi 
" 405aV^ 27a 4 Pl 2 „ + 5a$ 3a 2 + 2a 4 + 2 ' ^ J 

where di^ are two constants of integration with dimensions [di] = M~ 2 = \d<^ ■ It is interesting to note the 
presence of a new cosmological term in this Friedmann equation, which goes as a~ 5 , which would correspond to 
a matter term with equation of state parameter w = 2/3. 

If pa = 0,p m o = 0, i.e. a universe filled with radiation only, equation (82) has the following solution 

with [di] = [d 2 ] =M 2 . 

C. Case ci,C2 7^ 



In this case, one can divide equation (45) by c\ finding 



a A + c 3 p r0 + 3 a 4 / 
c 3 R — Qk 2 a 4 



where C3 = C2/C1 7^ 0. This implies that 



-3 2 

/flfl R = , 2 i r + -5— Ir H , (86) 

a(a^ + c 3 ) a z + c 3 
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where /i 3 , = /Xg/ci. 

Fricdmann equation Eq.(50) can be rewritten as 



f-^fu^-f R (R + 6 4)=-^- )+ ^ + ^. ( 87 ) 

a z + c 3 \ a z J a (a z + c 3 ) a 6 a 4 



By substituting (85) into (87), and solving for /, one finds 

12 Pl a 5 H (6k ~c 3 R) 



f 



a 4 (a 2 + c 3 ) [3 (12H 2 + R) a 4 + (30fc + c 3 R) a 2 + 18c 3 £;] 

p r o (12H 2 + R)a 4 + 2p m0 (c 3 R - 6k) a 3 + 3p r0 (c 3 R - 2k) a 2 + 6c 3 kp r0 
a 4 [3 (12H 2 + R)a 4 + (30fc + c 3 R) a 2 + I8c 3 k] 

which means that the Noether symmetry, combined with the dynamics, determines the form of /. In this case / is 
a function of o since both R and H are functions of a. We can still go further by using the same trick used in the 
previous section, i.e. considering / as an implicit function of a into the Noether condition (85). Since / = f(R(a)) 
one finds 

= d^_da_ dj_ _f_ 
lR dR dRda R'' { ' 

Plugging Eqs.(88) and (89) into (85), one finds a second order differential equation for H, as follows 

w 



a 4 (a 2 + c 3 ) (3a 2 + c 3 ) H 2 [18 plHa 3 + (a 2 + c 3 ) (4p r0 + Sap m0 )] 
x { -24c 3 (3a 2 + c 3 ) Mo H 4 a 5 - 24 (a 2 + c 3 ) 2 fc 2 ^ a 

-ff 2 [6 (3a 2 + c 3 ) 2 /ig i/' 2 a 4 + 24 (-3a 4 - 2c 3 a 2 + c 3 ) fc Mo 

+ (a 2 + c 3 ) 2 (45 p„ l0 a 3 + 72 p r0 a 2 + 21 c 3 p m0 a + 32 c 3 p r0 ) i?'] a 3 

- 6iJ 3 [(3a 2 + c 3 ) (15a 2 + 13c 3 ) /i 3 H' a 4 + 2c 3 (a 2 + c 3 f {2p r0 + ap m0 )] a 2 

- (a 2 + c 3 ) H [a 4 H' [12(C3 - 3a 2 ) k $ + (a 2 + c 3 ) (3a 2 + c 3 ) {4p r0 + 3ap m0 ) H'} 

- 4 (a 2 + c 3 ) fc (3p m0 a 3 + 6p r0 a 2 + 2c 3/ o r0 )] } . (90) 



This differential equation defines the dynamics of the Noether solutions for a generic f(R) model compatible with the 
Noether symmetry. This result is relevant since there is a free parameter c 3 , which together with the initial conditions 
for Hq and H' a , uniquely specify the dynamics. This non- linear ODE is still of second order in H (a) as the 0, O-Einstein 
equation for any f(R) theory. However, there is a huge improvement as this equation is independent of the explicit 
form f(R), having as the only unknown parameters two real numbers, c 3 and /zo, the Noether charge. This also says 
that for any value of the Noether charge there is a solution, the solution of (90). Therefore all the solutions of (90), 
as c 3 , po vary, represent the whole set of Noether-charged cosmological solutions of the f(R) theories. 



In this case equation (85) reduces to 



whereas / can be written as 



1. Vacuum and pure dust case 



, _. 3/(a 2 +c 3 ) 

fR -2(Rc 3 -6 K )> (91) 



, _ 2 (6k - Rc 3 ) ((6Hpl + p m p) a 2 + p m pc 3 ) 

•' ~ a (a 2 + c 3 ) (3 (12H 2 + R) a 4 + (30k + Rc 3 ) a 2 + 18kc 3 ) ' ^ ' 

The case p m o = 0, pa = admits no solutions, therefore, as before, we will only discuss the case po = 0, p m o ^ 0, for 
which we can recast / in the following form 

/ = 3 (12H 2 + R)a 4 + (30k + Rc 3 ) a 2 + 18kc 3 . (93) 
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Inserting this relation into (91) together with the definition of R one finds 



„ -4c 3 H 2 - a (15a 2 + 7c,) H' H - a 2 (3a 2 + c 3 ) H' A + 4k 
H" = i Sr^< uT J ( 94 ) 



a 2 (3a 2 + c 3 ) H 



whose general solution reads 



H-2 C3K 2K , 2dl _l 2c 3 rf 2 , „ , fQf ., 

2. Pure radiation case 

Once again, studying Eq. (90) to the case po = and p m o = 0, we find the following equation 

, 2 y,_ 18a 2 + 8c 3 , 2 y 12 c 3 H 2 2fc(6a 2 + 2c 3 ) 

1 j a (3a 2 + c 3 ) l j a 2 (3a 2 + c 3 ) a 4 (3a 2 + c 3 ) ' l J 

The general solution, when c3 > 0, for this ODE is 

2 3c 3 di 21d x 18di 5\Z3v^ d 2 9\/3d2 4k 2k 

a 4 c 3 or a' 3 a^/c 3 c 3 a z 

3c 3 d 2 arctan(^|) 27d 2 arctanf ^§f) 18d 2 arctanf^ff ) 
a 4 c 3 o z 

whereas, for c 3 < 0, one finds 

„ 2 3c 3 di 27rfi 18di hV?>Jc 3 'd 2 9V3d 2 4k 2k 

i? 2 = — -A + + -A- Xr 2 — + . - + — + — 



1 c 3 a" a° ffiy-C3 c 3 

3c 3 d 2 arctanh(^I§=) 27d 2 arctanh(^I§=) lSdsarctanhf^^ 



A3 



(98) 



Either expression for H(a) together with Eq. (88) and Eq. (60) form a solution for (9), and (44), and possess no = 
Nocther charge. 

3. Matter and Radiation case 
Let us restrict our study to the case p = 0, for which we can find analytical solutions. Eq.(90) reduces to 



. 2 \ " (45/9 m oa 3 + 72p r0 a 2 + 21c 3 p m0 a + 32c 3 p r0 ) , 2 w 

[H J = — - — [H J 

a (3a 2 + c 3 ) (4p r0 + 3ap m0 ) 

24c 3 (p m0 a + 2p r0 )g 2 8k{3p m0 a 3 + 6p r0 a 2 + 2c 3 p r0 ) 

a 2 (3a 2 + c 3 ) (4p r0 + 3ap m0 ) a 4 (3a 2 + c 3 ) (4p r0 + 3ap m0 ) 



(99) 



It is remarkable that this differential equation is linear in H 2 . This makes the problem of solving it much easier. In 
fact, analytical solutions for k = 0, ±1 can be achieved. Let us discuss them. 

1. k = 0. The solution of Eq.(99) is 

2 4did 2 c 9 3 /2 24did 2 cl' 2 po m d 2 c 5 3 /2 5/2 

ti = 1 5 V 6bdid 2 c a 



2V3p r0 arctan(^f)rf 2C 2 iQ Pr0 d 2C ?f 12x/3p r0 arctan(^f) d 2 c 3 



18 ^ 2 V^ + ISVWrctanf ^ ] d 2 , (100) 
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where d\ and d 2 are integration constants with dimensions, \d{\ — M 4 , and [o^] = M~ 2 . This is clearly a 
deviation from standard GR, because there is a 1/a term, which leads to an accelerated behavior if dominates. 
Furthermore there are terms, all involving p r0 , which include the arctangent of a, where c 3 is supposed to be 
positive. These terms have different behavior at low and high rcdshift. In fact since lim arctan(a) ~a at high 

redshifts, these terms behave as dust, 1/a and a respectively, and are subdominant with respect to the radiation. 
On the other hand, since lim arctan(a) ~ 7r/2 for large and positive a, these terms will behave as radiation, 

a— >oo 

curvature and cosmological constant respectively. It is also interesting to notice that in order to have a true 
dust matter component at late times, it has to be 

ttG 



10 p r0 d 2 C3 



3/2 



' pmO ■ 



(101) 



This means that p r o behaves as the source of matter component in this modified Fricdmann equation. A 
cosmological constant term is also present. It is determined by the integration constants of the Noether condition. 

As for the case c 3 < 0, the solution of Eq. (99) can be written as follows 



H 2 



4cM 2 (-c 3 ) 9 / 2 , 24d 1 d 2 (-c 3 y/ 2 P o m d 2 {-c 3 f/ 2 .. , , . , 5/2 
7 1 5 1 1 6Wia 2 {— c 3 ) 



2v3/9 r oarctanh( 



\/Za 



,.; d ^l , 10p r0 rf 2 (- C3 ) 3 / 2 



12\/3/5rO arc t an h( 



\/3a 



d 2 c 3 



18p r0 d 2 



(r 



18v3/9 r oarctanh 




(102) 



For this solution, as a pedagogical example, more detailed calculations and a link with scalar-tensor theories are 
given in the appendix. 



2. k / 0. The general solution is 
32k arctan f-^§=) 



H z 



n 3 

PrO 



9y/3a*fao^ 



16kp 2 q 
3a 2 p 2 re0 c 3 

9^2 PrO 

2ac| 



160K Pr0 
^a 3 p 3 m0 c 3 

_ 8kp% 
n 2 r 2 " T 

9%/3 arctan f-^|) d 2 p r0 



64k arctan 



(V3a\ 3 



8«PrO 

27oV 



3V3a 2 p 3 m0 4 /2 
a/3 arctan ( ^S ) d 2 p r0 



mX) 



2a 4 c 3 V2 
2k arctan (^ 



32kp 3 

3 «PmO C l 
5fi?2PrO 

2a 3 c| 



32k arctan (^f) p 3 

V^PmoCs 72 
3a/3 arctan (^=) d 2 p r0 



a^c 



7/2 



C3PrO 



4a/3k arctan f^S 



PrO 



1c 



9/2 



y/Stf 



PmO 



a Pm0\/C 3 



l0Kp r o QnprO 



3a 3 p m o ap m oC3 



6a/3k arctan ( ^|| J /Vo 
372 

PmOC 3 



2k 
ia 2 



K£3 
9a 4 



6di 9di di /9 m orf2 



a 2 c 3 



4a 4 c 2 



(103) 



Also in these cases we have interesting behaviors matching the main cosmological eras. The integration constants 
di j2 have dimensions respectively [d{\ = M 2 , and [d 2 ] — M~ 2 . The analysis, for both this and the previous case 
(k = 0), of the set of parameters {d\, d 2 ,c 3 } which can be bounded by observations will be done in a forthcoming 
paper. 

Eq. (99), for the case c 3 < 0, has solution 



H 2 = 



32Karctanh 



(7%) & 



64Karctanh 



( V3a \ 3 



9V3a 4 p 3 m0 , 



-c 3 

8np 2 r0 

3a 2 p 2 m0 c 3 27a 4 p 2 m0 



16np 



160kp 3 
27a 3 p 3 n0 c 3 3V3a 2 p 3 m0 (-c 3 ) 3 / 2 

V^arctanh (-)=) d 2 p 



32k P 3 



&Kp. 



rO 



n 2 c 2 
PmO c 3 



2a 4 (-c 3 ) 5 / 2 



9\/3arctanh ( -j== ) d 2 p r o 2Karctanh 



9d 2 p r0 

2ac 4 ~ ' 2(-c 3 ) 9 / 2 

10/sprO 6n Pr0 6V3Karctanh (-^l) p r0 



-y/3a 



32Karctanh (-7==) P r0 
3"P 3 mo4 + Vi Pm0 {~c 3 f/ 2 
M 2P ro , 3V3arctanh (^g|) rf 2 p r0 
2a 3 c 3 

-C 3 p r 



a 2 (-c 3 )V 2 
4A/3Karctanh (;-)=) pro 



V3a 4 p m o 



O?Pm0 



3a 3 p m0 ap m oC3 



Pmo(-C 3 ) 3/2 



2k 

3^ 



KC3 

9a 4 



6di 

a 2 C3 



9di 

^2~ 



dl 

„4 



-c?> 

Pmod 2 



(104) 
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It is worthy to note that once the free parameters are constrained by the data (the set of allowed parameters might 
be empty anyhow), one can select physically interesting f(R) models as in [21]. 

4- Non-linear case, po 7^ 

In this more general case, Eq.(90) cannot be written as a linear differential equation in H 2 , therefore it is not 
possible to achieve an analytical general solution. However, after fixing initial conditions for H and giving suitable 
values for the parameters, one can solve it numerically. These initial conditions fix, in turn, the f(R) model and the 
behavior of H{a). 

5. General non-linear case, c/0 and /io 7^ 

By using Eq. (43) inside Eq. (50) one finds the following expression for / 
CicR (12H 2 + R) a 5 cR (l2c 2 H 2 + 12ci« + c 2 R) a 3 



f 



(cia 2 + c 2 )A (cia 2 + c 2 )A 

2 (36ciK_ff/iQ — 6c2HR/j,q + 18cick 2 + 6c 2 Kp m o + 6c 2 c>ci? — ciC2p m oR) a 

(ci« 2 + c 2 )A 
2c 2 (— 18CK 2 — 6cip m QK + C 2 /9 m 0-R) PrO {\2c\H 2 a 4 + C\Ra 4 — 6cina 2 + 3c 2 i?a 2 + 6c 2 k) 



(cia 2 + c 2 )Aa A a 4 



(105) 



where 



A = 36ciiJ 2 a 4 + 3cii?a 4 + 30ci«a 2 + c 2 Ra 2 + 18c 2 k . (106) 



The Friedmann equation gives us the expression of / in terms of R(a), H{a) and a. Eq. (44), which can be rewritten 
here as 

f'(a) = 3 a (ci a 2 +c 2 ) /(a) - c (a 2 R(a) + 6k) ( Cl a 2 + c 2 ) p r0 

R'(a) 2a(c 2 R(a) - 6c lK ) 2a 4 (c 2 R(a) - 6cm) ' [ ' 

giving a dynamics for /, defines a second order differential equation for H, given by 

H" = [(cia 2 + c 2 ) HVy 1 H' 2 (l2c 2 CKa 7 + §c\p„ M a 7 + 54clfi 3 Ha 7 + 12c?p r0 a 6 + 2Acic 2 cKa 5 + 21c 2 c 2/ o„ l0 a 5 

+ 36cic 2 /i Ha + 28c 1 c 2 /? r oa + 12c 2 cKa + 15ciC 2 PmO« + 6c 2 p Ha + 20cic 2 p r aa + 3c 2 /9 TO oa + Ac^Pro) 

- [a (cia 2 + c 2 ) HT}' 1 H' (54c 2 cff 3 a 9 + 108cic 2 ci? 3 a 7 - 210c\^H 2 a 7 - 60c 2 c K ffa 7 - 45c 3 p rnQ Ha 7 

- 72c 3 p r0 Ha 6 + 36c 2 K^ 3 a 5 + 54c 2 ciJ 3 a 5 - 324cic 2 ^ 3 iJ 2 a 5 - 120cic 2 c«;iJa 5 - lllc 2 c 2/ 9 m0 ifa 5 

- 176c 2 c 2/ 9 r0 ifa 4 + 24cic 2 K^ 3 a 3 - 78c 2 /x 3 iJ 2 a 3 - 60c 2 cKi?a 3 - 87cic 2/ o m0 fla 3 - Yi6cic\p r QHa 2 

- \2c\np 3 a - 21clp m0 Ha - 32clp r0 H) - [a 4 ( Cl a 2 + c 2 ) i7 2 r] _1 4 (-18cic 2 ^ 3 ilV - 3c 2 c 2/ 9 m0 ff 3 a 7 
+ 18c 2 K^ 3 iJ 2 a 7 + 6c\cn 2 Ha 7 + Zc\np mQ Ha 7 - 6c 2 c 2 p r oH 3 a 6 + 6c\npraHa & - 6cjp, 3 H 4 a 5 - 6c 2 K 2 ^ 3 a 5 

- 6cic 2 p m o-ff 3 a 5 + 12ciC2Kp 3 H 2 a 5 + 12cic 2 cK 2 iJa 5 + 6c 2 C2np m oHa 5 — \2c\c\p r ^H 3 a + \Ac\c2np r oHa 

- 12c 1 c 2 K 2 p 3 a 3 - 3clp0mH 3 a 3 - 6c\K,p 3 H 2 a 3 + Qc\cK 2 Ha 3 + icic\np mn Ha 3 - 6c 3 2 p r vH 3 a 2 

+ lOcxcjnproHa 2 ~ Qcln 2 p 3 a + 2c 3 ,np ra H) , (108) 

where 

T = 18ciaff 2 a 7 + 18c 2 ci7 2 a 5 - 12cicna 5 - 9c 2 /9 m0 a 5 - bAcip 3 Ha b - I2c\p r0 a 4 - 12c 2 cna 3 - 12cic 2/ o m0 a 3 

— lSc2pf Ha — 16cic 2 p r oa — 3c 2 p m o a — 4c 2/ o r o . (109) 

It is evident that a more detailed (numerical) study, pursued elsewhere, of this differential equation is necessary in 
order to study the dynamics of these solutions. 
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D. Non-Noether solutions 

In general it is not possible to find a solution of the Friedmann equations which is also a Noether symmetry since, 
in principle, such symmetries do not exist for any f(R) theory. In general, a solutions of the cosmological equations 
is not a solution compatible with the condition L-^C = 0. This is a peculiar situation which holds only if conserved 
quantities (Noether's charges) are intrinsically present in the structure of the theory (in our case, the form of ,f(R))- 
For example, imposing a power law solution, a ex t p , defines a function of R — R(a), which can be put in the Noether 
symmetry equations, in order to find / = f(R(a)). Finally one can substitute the expressions for /(a), R(a), and H 
in the Friedmann equations. In doing this, it is easy to show that, for k — 0, there are no simple power-law solutions 
compatible with a Noether charge. 

The method discussed above allows to discriminate theories which admit or not cosmological solutions compatible 
with a Noether charge. 

It is also clear that power-law solutions do exist in general for f(R) models, but they can be found using different 
methods [24]. Assuming, in general, a power-law H(a), one finds R as a function of a, and then, in principle, 
/ = f(R(a)). It is therefore possible to write the Einstein equation as a second order differential equation for / as 
a function of a, whereas all other quantities (H and R) are given functions of a. The same argument holds for the 
redshift z [21]. 

For example, let us rewrite the Friedmann equation (8) as 

f-6f lm RH-6f R H 2 -f R (R+^)=^ + ^, (110) 

\ a z ) a J a 4 

and let us consider H = H(a) and R = R{a) as given functions of a, being, as above, 

- - k 

R= -12H 2 -6aHH' -6^. (Ill) 



a 2 



The Friedmann equation can be written as 



1 R" 1 /- 6k 

R+^T 



a R' 6a H 2 



f-^L-f = - Pm0a + P /° Ri. (112) 

J 6aH 2J 6a 5 H 2 V ' 



This is a second order linear equation in /, whose general solutions depends on two parameters, /o and /q. Specifically, 
being the equation linear, the general solution is the linear combination of two solutions of the homogeneous ODE 
plus a particular solution. It is then clear that more than one f(R) model can have the same behavior for H(a), i.e. 
more theories share the same cosmological evolution. This situation is due to the fact that one has a fourth-order 
gravity theory. The singular points of this differential equation are those for which either H or dR/da vanishes. 
Starting from these considerations, interesting classes of solutions can be found out. 

1. Radiation solutions 
Let us seek for all the f(R) models which have the particular solution a = yt/to, which means 



1 H , -6k 

_ ~ = —^ , so that R = 5 

2 to a or a z 



H = T7TT2 = TT ' so that R = ^ , (113) 



where Hq = (2 to) 1 . We have three interesting cases. 

1. For k — 0, we have R = 0, leading to the Friedmann equation 

f(Q)-6f R (0)H 2 = ^ + ^, (114) 

a 6 a 

which, if p m o y^ 0, cannot be solved for H ~ aT 2 since /(0) and /'(0) cannot be functions of a, but only 
constants. If p m o = 0, standard GR is of course recovered. 

2. For the case k — — 1 we have the following differential equation for /, 

f ,i , 4 ,,, 2k f 2 k (p r0 + ap mQ ) 

1 + a f+ Hl f = Hja-i ' (U5) 



whose general solution can be written as 



6k 

-ft — 7T 

a z 



f = 



/ a\f — k j „ n „( a-J —2k \ rrl / a\/ — k j „■ / a-J —2k \ tj2 4 /o / a\f — k j „ , _ / a\/ — 2k \ tj 



^2a 7 /2 K% /^F 



^2a 7 /2 K ^f 



7,5/2,/---; 



4/77 / a-J — k j r ■ / a-v/— 2k \ ft / /-1 • / \/2 a v/ — k \ ,* / a-J — 2k \ 

V2 A /^-d 2 sm^^— jff„ pm0 ProV^Ci^^-j*— J sin^-^J 



PrOV — kcos 



a\f—2K\ n: / a\f —Ik 



Ho 



Si 



Ho 



p r K cos 



Ho 



j2a 3 H Q 



Ci 



Ho 



\/2a 3 #o 

Ho—) Sl { Ho— ) 



a 2^2 



where the Sinlntegral and Coslntcgral functions, Si and Ci respectively, are defined as 



Si(,. )= / ^±dt Ci(a;) = - 



cos(t) 



dt. 



The integration constants dip have dimensions [di] = [d 2 ] = Af 4 . 
3. Along the same lines, the case k — 1 has the following solution 






/ = 



^E dl cosh ^ ff 2 ^ dl sinh ±ff fl* tffy^dx cosh ^ H Q 



J2a 7 / 2 K J^ 



j2aV 2 nJ^ 



^^d l8 mh(^)Fo pm0 p,o^Chi(^)sinh(±g 



2kq, 



I 5 / 2 , 



/2/sa \ cv,i / V2ko \ - _ _- u / \f2na \ /"iu; / V^kcl 



p r0 ^cosh( ^f j Shif ^f Pr0 «coshf ^f Chif ^ 



j2a 3 H a 



a 2 H 2 



Pr0KS inh(^jShi(^ 



where the hyperbolic Sinlntegral and Coslntegral, Shi and Chi respectively, are defined as 

Shi(x) : 



x sinh(i) , „,., . , - > r cosh(t) - 1 , 
; dt Chi(x) = 7b,m + ln(a;) + / — dt . 



t 



t 
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(116) 



(117) 



(118) 



(119) 



(120) 



(121) 



and 7b m ~ 0.577 is the Euler-Mascheroni constant. Both d\ and d 2 are integration constants which dimensions 

M 4 . 



2. Matter solutions 

In this case, we search for f(R) models which have a dust-matter behavior, that is a — (t/to) 2 ' 3 , 

2 H a , - 2(2/4 + 9 ka) 



H 



and R 



3t a 3 / 2 a 3/2 ' 3a 3 

where Hq = 2/ (3 to). For the case k = 0, we find the explicit analytic solution 

4 



R = -- 



3 1(, a 



2„3 ' 



./(a) = a -( 7 +^)/ 4 (di a^ /2 + da) + 



PmO a — 6p r o 

ifo 1 



(122) 

(123) 
(124) 
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This is a 2-parameters family of solutions, depending on the two integration constants di,2 both with dimensions M 4 . 
The Einstein- Hilbert case f(R) = R belongs to this family, when di, d 2 , and p r o all vanish. 



3. Exponential solutions 
In this case, we look for the behavior 

H = Hq = constant , which is 

As above, we have three cases depending on k. 



r = -vim- 



()k 



1. k = 0. Both H and R are constants, and R = Rq = — 12Hq. The Friedmann equation is 



f(Ro) — t fii(Ro) Ro 



PmO PrO 



and it has solutions only for p m Q = p r Q = being Ro a constant (see also [53]). 
2. k = 1. In this case, H is still a constant but R is not. One finds 

6 K 



R = 


-- l2H «-^ 


I -- 


= d\ cosh — — + d 2 smh — — 
\H aJ yH Q aJ 




6p r o H 3p m o H 6p r o H p r Q p m o 

9 9 4 3 



3. fc 



-1. The solution is 



6 k 



R = 


- - UH 0-^2> 


f = 


= «i cos — + d 2 sin — 

V H Q a ) \ H a ) 

6p r0 H Q 'ipmQ H 6p r0 H p r0 p m0 

9 9 4 ! 3 



(125) 



(126) 



(127) 



(128) 



(129) 



(130) 



4. ACDM solutions 

Let us now look for f(R) models which are compatible with the ACDM being solutions of Friedmann equations. 
This analysis could be extremely important to compare the f(R) approach with observations (see also [47]). One 
defines 



H 2 = Hi 



n 



mO 



fl r Q 



1 — "mil — " 



a a 

The differential equation to solve is therefore the following 

6fl m ()H 4(fi m o + ^K> 



M) 



/" 



Da 4 



7\ IruO® — 8'2 r Q 



3^ m0 7?Q 4- Aak — (fimO + ^ro — l)a + ^mOfl + ^ r o 



2« 



3fi m0 #n + 4afc 



/ 



2a [-(fi m o + Qro - 1)« 4 + ^moa + ^rf>] -ffo 

(3^ m0 gg + 4ak) (p r0 + ap m0 ) 

2a 5 [—(fl m o + ^ro - 1)« 4 + f^moa + ^ro] Hq 



(131) 



(132) 



The general integral can be numerically achieved by giving suitable initial conditions for / , /q. This analysis will be 
pursued in a forthcoming paper. 
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VI. DISCUSSION AND CONCLUSIONS 

In this paper, we have discussed a general method to find out exact/analytical cosmological solutions in f(R) 
gravity. The approach is based on the search for Noether symmetries which allow to reduce the dynamics and, in 
principle, to solve more easily the equations of motion. Besides, due to the fact that such symmetries are always 
related to conserved quantities, such a method can be seen as a physically motivated criterion. 

The main point is that the existence of the symmetry allows to fix the form of f(R) models assumed in a point-like 
cosmological action where the FLRW metric is imposed. It is worth noticing that, starting from a point-like FLRW 
Lagrangian, and then deriving the Euler-Lagrange equations of motion, leads exactly to the same equations obtained 
by imposing the FLRW metric in the Einstein field equations. This circumstance allows to search "directly" the 
Noether symmetries in the point-like Lagrangian and then to plug the related conserved quantities into the equations 
of motion. As a result i) the form of the f{R) is fixed directly by the symmetry existence conditions and ii) the 
dynamical system is reduced since some of its variables (at least one) is cyclic. 

The method is useful not only in a cosmological context but it works, in principle, every time a canonical, point-like 
Lagrangian is achieved (in [54], it has been used to find out spherically symmetric solutions in f(R) gravity). 

In this paper, we have considered a generic f(R) theory where standard fluid matter (dust and radiation) is present. 
The Noether conditions for symmetry select forms of f(R) depending on a set of cosmological parameters such as 
{Pro, pmOi k, Ho} and the effective gravitational coupling. Such a dependence can be easily translated into the more 
suitable set of observational parameters {Sl r o, U ra o, ^fc, Ho} and then matched with data. This situation has a twofold 
relevance: from one side, it could contribute to remove the well known problem of degeneracy (several dark energy 
models fit the same data and, essentially, reproduce the ACDM model); from the other side, being the search for 
Noether symmetries a relevant approach to find out conserved quantities in physics, this could be an interesting 
method to select models motivated at a fundamental level. It is worth noticing that the Noether constant of motion, 
which we have found, has the dimensions of a mass and is directly related to the various sources present into dynamics. 
In some sense, the Noether constant "determines" the bulk of the various sources as p m o 1 p r o and the effective p\ 
and then could greatly contribute to solve the dark energy and dark matter puzzles. In a forthcoming paper, we will 
directly compare the solutions which we have presented here with observational data. 

The "non- Noether solutions" deserve a final remark. In this case, we do not ask for a Noether symmetry but, 
finding these solutions, can be related to the previous general method. We have shown that the standard cosmological 
behaviors of the usual Einstein-Fricdmann cosmology can be achieved also in generic f(R) models, assuming that 
the cosmological quantities H and R depend on the scale factor a. As result, we find out general f(R(a)) where the 
standard solutions of the linear f(R) = R case are easily recovered. 

Acknowledgment. We want to thank prof. Ringeval and prof. Fabri for useful discussions and comments. ADF is 
supported partly by STFC, UK and partly by the Belgian Federal Office for Scientific, Technical and Cultural Affairs 
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Appendix A: Solutions and link with scalar-tensor theories 

We will explicitly show, as an example, that equation (102) is indeed a Noether solution (with k — 0, flat space, 
and no = 0, zero Noether charge). First, from H(a), given by 



H 2 = - 



4did 2 (-c 3 ) 9 / 2 24d 1 d 2 (-c 3 ) 7 / 2 p 0m d 2 (-c 3 ) 5 / 2 /2 



<r <r 



2 v / 3/5roarctanh( : ^jd2c| I0p r od 2 (-c 3 ) 3 / 2 12V3/9 r oarctanh( -^=) d 2 c 3 







(Al) 
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(A2) 



we can calculate the expression for R(a) as follows 
R = -12 H 2 - 6 aHH' 

144rf 1 rf 2 (-c 3 ) 7/2 , , ^ , , , 5/2 48rf 2/ yo(-c 3 ) 3 / 2 72V3d 2 p r0 axctanh(^) c 3 

H ■ — 216v3(i2/O r o arc tanh — == 

Since we know both _ff and R, now, by using Eq. (88), we can find f(a) as follows 

8d lC § 24rfiC3 3p r0 4V3p r0 arctanh(^) i 2 p r0 12^ r0 arctanh(^) 
a a a* a^y/^c^ a z c 3 a(— c 3 ) J ' z 

These expressions for /, _R, H fulfill equation (9). The system has also a constant of motion /uo = given by equation 
(48), as the Lagrangian possesses a Noether symmetry. 

We will discuss how to link this solution (extending this procedure to the other solutions is straightforward) to the 
scalar-tensor picture, by finding the potential for the scalar non-minimally coupled with gravity. In fact, starting from 
the action 

S = J d A x^f{R) + S m , (A4) 

one can rewrite it (at least at the classical level) in the following form 

S = I d A xsP^{U R - V(cp)} + S m , (A5) 



where V = cp f v — f(<p), and f v — df/dtp. The classical equation of motion for <p leads to <p — R. One can make a 
field redefinition to bring the action in the form 

<£xyj—g[-x R - V(x)} + S m , (A6) 

where x = ~f<p- 

In this case we can use our solutions in order to find V(x), the only unknown in the theory. One can do it as follows 

X = -U = ~fn = ~ (A7) 

V = ?U-f = Rf R -f = R^-R, (A8) 

where these relations are correct on shell, i.e. for the solutions of the equations of motion. Using equations (Al), 
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(A2), and (A3), one can write down explicitly the potential, at least for this case, as follows 
V(x) = 3456di^ x 3(_ C3 )i3/ 2 _ 103684^0 X 4 cf 



1728^3^/9,-0 x 3 arctanh 



\/3 6d 2 x(-c 3 ) 5/2 + V-36^x 2 ^-c 3 



+ 17284p r0 X 3 \f-36dlx 2 4-c 3 (-c 3 ) 7 / 2 - 288d lC ? 2 x(-c 3 ) 5 / 2 + 432d 2 p r0 X 2 c; 



288v / 3dl/9 r0 X 2 \/-36d| x 2 c| - c 3 arctanh 



V3 6d 2 x(-c 3 ) 5/2 + V-^4x 2 4-c 3 



(~c 3 ) 3 / 2 



144V3<i 2 p r o X arctanh 



V3 (6d 2 x(-c 3 ) 5/2 + V-36^x 2 ^-c 3 



16d n /-36d|x 2 c|-c 3 



96d 2 PrO X-\/-36rf2 X 2 c 3 - c 3 9p r0 



- - 576cM 2 x 2 V-36^ x 2 c| - c^ 



+ 8VT5p r0 (-c,)-^ 2 ^-36d 2 2 x 2 4 - ca 



arctanh 



V3 (6d 2 x(-c 3 ) 5/2 + V-36^ x 2 cl-c 3 



(A9) 



In order to study the evolution of the background, whether or not it leads to a viable dynamics for the universe, it is 
already sufficient to check if the Hubble parameter given by (Al) can fit the data, from Big Bang Nucleosynthesis up 
to Dark Energy domination. 
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